################################################################################################################
################################################################################################################
#### Levels of Linkage: Across-Agreement v. Within-Agreement Explanations of Consensus Formation Among States   
####   Heather McKibben & Shaina Western #######################################################################
################################################################################################################
#################################################################################################################

library(foreign)
library(MCMCpack)
extractCodaDraws <- function(x, vars=varnames(x)) {
  # Checks
  if( class(x) != "mcmc.list" ) {
    stop("extractCodaDraws requires an object of class 'mcmc.list'")
  }
  for(vn in vars) {
    if( !(vn %in% varnames(x)) ) {
      stop("Each varible selected must exist in the 'mcmc.list'")
    }
  }
  # extract draws
  out <- sapply(vars, function(vn) unlist(lapply(x, function(tc) tc[,vn])))
  rownames(out) <- NULL
  return( as.data.frame(out) )
}


### Load Data ###
DataSet <- read.csv("/Dataset for EU Consensus behavior McKibben-Western 2010-10-17-5.csv")
DataSet$LogPolicyLink <- log(DataSet$Policy.Macro.Linkage)
DataSet$LogTemporalLink <- log(DataSet$Temporal.Macro.Linkage)
DataSet$Consensus <- c(rep(NA, 70))
for(i in 1:length(DataSet$Consensus)){
	DataSet$Consensus[i] <- ifelse(DataSet$QMV[i]==0, 1, 0)
}

DataSet$LogCouncilLink <- log(DataSet$p0 + 1)

### Alternative measurement of missing 1

Model.ods1a <- MCMClogit(Consensus2 ~ Microlink_ods1 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=2590)

Model.ods1b <- MCMClogit(Consensus2 ~ Microlink_ods1 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=679)

Model.ods1c <- MCMClogit(Consensus2 ~ Microlink_ods1 + LogPolicyLink +  LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=9145)

Model.ods.1 <- mcmc.list(Model.ods1a, Model.ods1b, Model.ods1c)



Model.ods2a <- MCMClogit(Consensus2 ~ Microlink_ods2 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=2590)
	
Model.ods2a <- MCMClogit(Consensus2 ~ Microlink_ods2 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=679)

Model.ods2a <- MCMClogit(Consensus2 ~ Microlink_ods2 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=9145)

Model.ods2 <- mcmc.list(Model.ods2a, Model.ods2b, Model.ods2c)


Model.ods3a <- MCMClogit(Consensus2 ~ Microlink_ods3 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=2590)

Model.ods3b <- MCMClogit(Consensus2 ~ Microlink_ods3 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=679)

Model.ods3c <- MCMClogit(Consensus2 ~ Microlink_ods3 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=9145)

Model.ods3 <- mcmc.list(Model.ods3a, Model.ods3b, Model.ods3c)

Model.ods4a <- MCMClogit(Consensus2 ~ Microlink_ods4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=2590)

Model.ods4b <- MCMClogit(Consensus2 ~ Microlink_ods4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=679)

Model.ods4c<- MCMClogit(Consensus2 ~ Microlink_ods4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=9145)

Model.ods4 <- mcmc.list(Model.ods4a, Model.ods4b, Model.ods4c)

Model.ods5a <- MCMClogit(Consensus2 ~ Microlink_ods4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=2590)

Model.ods5b <- MCMClogit(Consensus2 ~ Microlink_ods4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=679)

Model.ods5c <- MCMClogit(Consensus2 ~ Microlink_ods4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=9145)

Model.ods5 <- mcmc.list(Model.ods5a, Model.ods5b, Model.ods5b)

summary(Model.ods1)
summary(Model.ods2)
summary(Model.ods3)
summary(Model.ods4)
summary(Model.ods4)

DrawModel.ods1 <- extractCodaDraws(Model.ods1)
DrawModel.ods2 <- extractCodaDraws(Model.ods2)
DrawModel.ods3 <- extractCodaDraws(Model.ods3)
DrawModel.ods4 <- extractCodaDraws(Model.ods4)
DrawModel.ods5 <- extractCodaDraws(Model.ods4)


mean(DrawModel.ods1$Microlink_ods1>0)
mean(DrawModel.ods2$Microlink_ods2>0)
mean(DrawModel.ods3$Microlink_ods3>0)
mean(DrawModel.ods4$Microlink_ods4>0)
mean(DrawModel.ods5$Microlink_ods5>0)

mean(DrawModel.ods1$LogPolicyLink>0)
mean(DrawModel.ods2$LogPolicyLink>0)
mean(DrawModel.ods3$LogPolicyLink>0)
mean(DrawModel.ods4$LogPolicyLink>0)
mean(DrawModel.ods5$LogPolicyLink>0)

mean(DrawModel.ods1$LogTemporalLink>0)
mean(DrawModel.ods2$LogTemporalLink>0)
mean(DrawModel.ods3$LogTemporalLink>0)
mean(DrawModel.ods4$LogTemporalLink>0)
mean(DrawModel.ods5$LogTemporalLink>0)

mean(DrawModel.ods1$LogCouncilLink>0)
mean(DrawModel.ods1$LogCouncilLink>0)
mean(DrawModel.ods1$LogCouncilLink>0)
mean(DrawModel.ods1$LogCouncilLink>0)
mean(DrawModel.ods1$LogCouncilLink>0)


### Alternative measurement of missing 2

Model.odb1a <- MCMClogit(Consensus2 ~ Microlink_odb1 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=2590)

Model.odb1b <- MCMClogit(Consensus2 ~ Microlink_odb1 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=679)

Model.odb1c <- MCMClogit(Consensus2 ~ Microlink_odb1 + LogPolicyLink +  LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=9145)

Model.odb.1 <- mcmc.list(Model.odb1a, Model.odb1b, Model.odb1c)



Model.odb2a <- MCMClogit(Consensus2 ~ Microlink_odb2 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=2590)

Model.odb2a <- MCMClogit(Consensus2 ~ Microlink_odb2 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=679)

Model.odb2a <- MCMClogit(Consensus2 ~ Microlink_odb2 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=9145)

Model.odb2 <- mcmc.list(Model.odb2a, Model.odb2b, Model.odb2c)


Model.odb3a <- MCMClogit(Consensus2 ~ Microlink_odb3 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=2590)

Model.odb3b <- MCMClogit(Consensus2 ~ Microlink_odb3 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=679)

Model.odb3c <- MCMClogit(Consensus2 ~ Microlink_odb3 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=9145)

Model.odb3 <- mcmc.list(Model.odb3a, Model.odb3b, Model.odb3c)

Model.odb4a <- MCMClogit(Consensus2 ~ Microlink_odb4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=2590)

Model.odb4b <- MCMClogit(Consensus2 ~ Microlink_odb4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=679)

Model.odb4c<- MCMClogit(Consensus2 ~ Microlink_odb4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=9145)

Model.odb4 <- mcmc.list(Model.odb4a, Model.odb4b, Model.odb4c)

Model.odb5a <- MCMClogit(Consensus2 ~ Microlink_odb4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=2590)

Model.odb5b <- MCMClogit(Consensus2 ~ Microlink_odb4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=679)

Model.odb5c <- MCMClogit(Consensus2 ~ Microlink_odb4 + LogPolicyLink + LogTemporalLink + LogCouncilLink + 
	Codecision + Consensus + Regulation, data=DataSet, burnin=10000000, mcmc=20000000, 
	thin=10000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=9145)

Model.odb5 <- mcmc.list(Model.odb5a, Model.odb5b, Model.odb5b)

summary(Model.odb1)
summary(Model.odb2)
summary(Model.odb3)
summary(Model.odb4)
summary(Model.odb4)

DrawModel.odb1 <- extractCodaDraws(Model.odb1)
DrawModel.odb2 <- extractCodaDraws(Model.odb2)
DrawModel.odb3 <- extractCodaDraws(Model.odb3)
DrawModel.odb4 <- extractCodaDraws(Model.odb4)
DrawModel.odb5 <- extractCodaDraws(Model.odb4)


mean(DrawModel.odb1$Microlink_odb1>0)
mean(DrawModel.odb2$Microlink_odb2>0)
mean(DrawModel.odb3$Microlink_odb3>0)
mean(DrawModel.odb4$Microlink_odb4>0)
mean(DrawModel.odb5$Microlink_odb5>0)

mean(DrawModel.odb1$LogPolicyLink>0)
mean(DrawModel.odb2$LogPolicyLink>0)
mean(DrawModel.odb3$LogPolicyLink>0)
mean(DrawModel.odb4$LogPolicyLink>0)
mean(DrawModel.odb5$LogPolicyLink>0)

mean(DrawModel.odb1$LogTemporalLink>0)
mean(DrawModel.odb2$LogTemporalLink>0)
mean(DrawModel.odb3$LogTemporalLink>0)
mean(DrawModel.odb4$LogTemporalLink>0)
mean(DrawModel.odb5$LogTemporalLink>0)

mean(DrawModel.odb1$LogCouncilLink>0)
mean(DrawModel.odb1$LogCouncilLink>0)
mean(DrawModel.odb1$LogCouncilLink>0)
mean(DrawModel.odb1$LogCouncilLink>0)
mean(DrawModel.od1$LogCouncilLink>0)


##########################################################################################################
##########################################################################################################
### Frequentist Model
##########################################################################################################
##########################################################################################################

FrequentistModel <- glm(Consensus2 ~ Microlink_ODMB1 + LogPolicyLink + LogTempralLink + LogCouncilLink + 
	Codecision + QMVInv + Regulation, data=DataSet, family=binomial(link="logit"))
FrequentistModel <- glm(Consensus2 ~ Microlink_ODMB2 + LogPolicyLink + LogTempralLink + LogCouncilLink + 
	Codecision + QMVInv + Regulation, data=DataSet, family=binomial(link="logit"))
FrequentistModel <- glm(Consensus2 ~ Microlink_ODMB3 + LogPolicyLink + LogTempralLink + LogCouncilLink + 
	Codecision + QMVInv + Regulation, data=DataSet, family=binomial(link="logit"))	
FrequentistModel <- glm(Consensus2 ~ Microlink_ODMB4 + LogPolicyLink + LogTempralLink + LogCouncilLink + 
	Codecision + QMVInv + Regulation, data=DataSet, family=binomial(link="logit"))
FrequentistModel <- glm(Consensus2 ~ Microlink_ODMB5 + LogPolicyLink + LogTempralLink + LogCouncilLink + 
	Codecision + QMVInv + Regulation, data=DataSet, family=binomial(link="logit"))
	
summary(freq1)
summary(freq2)
summary(freq3)
summary(freq4)
summary(freq5)


##########################################################################################################
##########################################################################################################
#Get the data
DataSetNoMiss <- read.csv("NoMissing.csv")
DataSet$LogPolicyLink <- log(DataSet$Policy.Macro.Linkage)
DataSet$LogTempralLink <- log(DataSet$Temporal.Macro.Linkage)
DataSet$Consensus <- c(rep(NA, 70))
for(i in 1:length(DataSet$Consensus)){
	DataSet$Consensus[i] <- ifelse(DataSet$QMV[i]==0, 1, 0)
}

DataSet$LogCouncilLink <- log(DataSet$p0 + 1)

DataSet$p0 <- DataSet$p0 + 1
DataSet$logp0 <- log(DataSet$p0)


library(MCMCpack)

NoMissinga <- MCMClogit(Consensus2 ~ MicrolinkageNoMiss + LogPolicyLink + LogTempralLink + logp0 + Codecision + QMVInv + Regulation, data= DataSetNoMiss, burnin=10000000, mcmc=5000000, thin=2000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=2590) 

NoMissingb <- MCMClogit(Consensus2 ~ MicrolinkageNoMiss + LogPolicyLink + LogTempralLink + logp0 + Codecision + QMVInv + Regulation, data= DataSetNoMiss, burnin=10000000, mcmc=5000000, thin=2000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=679)

NoMissingc <- MCMClogit(Consensus2 ~ MicrolinkageNoMiss + LogPolicyLink + LogTempralLink + logp0 + Codecision + QMVInv + Regulation, data=DataSetNoMiss, burnin=10000000, mcmc=5000000, thin=2000, tune=.7, B0=.00001, b0=0, marginal.likelihood="Laplace", seed=9145)

NoMissing <- mcmc.list(NoMissinga, NoMissingb, NoMissingv)
summary(NoMissing)
DrawNoMissing <- extractCodaDraws(NoMissing)
mean(DrawNoMissing $MicrolinkageNoMiss>0)
